Basin attractors for various methods

M. Scott et al. / Applied Mathematics and Computation 218 (2011) 2584–2599 2585 (1) How does the basin of attraction differ for algorithms with the same order of convergence. (2) How does the basin of attraction differ for algorithms with different order of convergence. (3) Can the differences be used to compare various algorithms? In this paper we will discuss some qualitative issues using the basin of attraction as a criterion for comparison. To this end, we shall recall some preliminaries, see for example Milnor [9] and Amat et al. [1]. Let R : b C ! b C be a rational map on the Riemann sphere. Definition. For z 2 b C we define its orbit as the set orbðzÞ 1⁄4 fz;RðzÞ;RðzÞ; . . . ;RðzÞ; . . .g: Definition. A point z0 is a fixed point of R if R(z0) = z0. Definition. A periodic point z0 of period m is such that R(z0) = z0 where m is the smallest such integer. Remark 1. If z0 is periodic of period m then it is a fixed point for R. We classify the fixed points of a map based on the magnitude of the derivative. Definition. A point z0 is called attracting if jR(z0)j < 1, repelling if jR(z0)j > 1, and neutral if jR(z0)j = 1. If the derivative is also zero then the point is called super-attracting. Definition. The Julia set of a nonlinear map R(z), denoted J(R), is the closure of the set of its repelling periodic points. The complement of J(R) is the Fatou set FðRÞ: By its definition, J(R) is a closed subset of b C. A point z0 belongs to the Julia set if and only if dynamics in a neighborhood of z0 displays sensitive dependence on the initial conditions, so that nearby initial conditions lead to wildly different behavior after a number of iterations. As a simple example, consider the map R(z) = z on b C. The entire open disk is contained in FðRÞ; since successive iterates on any compact subset converge uniformly to zero. Similarly the exterior is contained in FðRÞ: On the other hand if z0 is on the unit circle then in any neighborhood of z0 any limit of the iterates would necessarily have a jump discontinuity as we cross the unit circle. Therefore J(R) is the unit circle. Such smooth Julia sets are exceptional. Invariance Lemma [9]: The Julia set J(R) of a holomorphic map R : b C ! b C is fully invariant under R. That is, z belongs to J if and only if R(z) belongs to J. Iteration Lemma: For any k > 0, the Julia set J(R) of the k-fold iterate coincides with J(R). Definition. If O is an attracting periodic orbit of period m, we define the basin of attraction to be the open set A 2 b C consisting of all points z 2 b C for which the successive iterates R(z), R(z), . . . converge towards some point of O. Lemma 1. Every attracting periodic orbit is contained in the Fatou set of R. In fact the entire basin of attraction A for an attracting periodic orbit is contained in the Fatou set. However, every repelling periodic orbit is contained in the Julia set. The idea of basin of attraction of some root-finding methods was introduced by Stewart [17]. He compared Newton’s method to the third order methods given by Halley [4], Popovski [5] and Laguerre [6]. In an ideal case, if a function has n distinct zeros, then the plane is divided to n basins. For example, if we have the polynomial z 1, then the roots are z = 1 and z 1⁄4 1 ffiffi 3 p i 2 , see Fig. 1. Ideally the basins boundaries are straight lines. Actually, depending on the numerical method, we find the basin boundaries are much more complex, see examples later. Our study considers ten methods of various orders, two of which were considered by Stewart [17]. We include optimal methods of order p = 2,4,8,16. Note that a method of order p = 2 is optimal (see [8]) in the sense that it requires n + 1 function-(and derivative-) evaluations per cycle. The methods we consider here with their order of convergence are: (1) Newton’s optimal method (p = 2). (2) Halley’s method (p = 3). (3) King’s family of optimal methods (p = 4). (4) Kung–Traub’s optimal method (p = 4). (5) Murakami’s method (p = 5). (6) Neta’s family of methods (p = 6). (7) Chun–Neta’s method (p = 6). (8) Neta–Johnson’s method (p = 8). (9) Neta–Petkovic’s optimal method (p = 8). (10) Neta’s family of optimal methods (p = 16). –0.8 –0.6 –0.4 –0.2 0 0.2 0.4 0.6 0.8 –0.4 –0.2 0.2 0.4 0.6 0.8 1 Fig. 1. Location of the roots of z 1. 2586 M. Scott et al. / Applied Mathematics and Computation 218 (2011) 2584–2599 The reason why we introduced more than one optimal fourth order method and more than one sixth order method will be clarified later. Newton’s optimal method (see e.g. Conte and deBoor [3]) is of second order for simple roots and given by xnþ1 1⁄4 xn un; ð1Þ where un 1⁄4 fn f 0 n ð2Þ and fn = f(xn) and similarly for the derivative. Halley’s method [4] is of third order and given by xnþ1 1⁄4 xn un 1 f 00 n 2f 0 n un : ð3Þ King’s fourth order optimal family of methods [7] is given by yn 1⁄4 xn un; xnþ1 1⁄4 yn f ðynÞ f 0 n fn þ bf ðynÞ fn þ ðb 2Þf ðynÞ : ð4Þ For the case b = 0 the method is actually due to Ostrowski [16]. Another optimal fourth order method is due to Kung and Traub [8] given by yn 1⁄4 xn un; xnþ1 1⁄4 yn f ðynÞ f 0 n 1 1⁄21 f ðynÞ=fn 2 : ð5Þ Murakami’s fifth order method [10] is given by xnþ1 1⁄4 xn a1un a2w2ðxnÞ a3w3ðxnÞ wðxnÞ; ð6Þ where un is given by (2) and w2ðxnÞ 1⁄4 fn f 0ðxn unÞ ; w3ðxnÞ 1⁄4 fn f 0ðxn þ bun þ cw2ðxnÞÞ ; wðxnÞ 1⁄4 fn b1f 0 n þ b2f 0ðxn unÞ : ð7Þ To get fifth order, Murakami suggested several possibilities and we picked the following M. Scott et al. / Applied Mathematics and Computation 218 (2011) 2584–2599 2587 c 1⁄4 0; a1 1⁄4 :3; a2 1⁄4 :5; a3 1⁄4 2 3 ; b1 1⁄4 15 32 ; b2 1⁄4 75 32 ; b 1⁄4 1 2 : ð8Þ Neta’s sixth order family of methods [11] is given by yn 1⁄4 xn un; zn 1⁄4 yn f ðynÞ f 0 n fn þ bf ðynÞ fn þ ðb 2Þf ðynÞ ; xnþ1 1⁄4 zn f ðznÞ f 0 n fn f ðynÞ fn 3f ðynÞ : ð9Þ Note that the first two substeps are King’s method. Several choices for the parameter b were discussed. Chun and Neta [2] show that b 1⁄4 2 is best. Another sixth order method due to Chun and Neta [2] is based on Kung and Traub scheme [8], yn 1⁄4 xn un; zn 1⁄4 yn f ðynÞ f 0 n 1 1⁄21 f ðynÞ=fn 2 ; xnþ1 1⁄4 zn f ðznÞ f 0 n 1 1⁄21 f ðynÞ=fn f ðznÞ=fn 2 : ð10Þ Neta and Johnson [12] have developed an eighth order method based on Jarratt’s method [13] yn 1⁄4 xn un; zn 1⁄4 xn fn 1 6 f 0 n þ 6 f ðynÞ þ 3 f ðgnÞ ; gn 1⁄4 xn 1 8 un 3 8 fn f ðynÞ ; xnþ1 1⁄4 zn f ðznÞ f 0 n f 0 n þ f ðynÞ þ a2f ðgnÞ ð 1 a2Þf 0 n þ ð3þ a2Þf ðynÞ þ a2f ðgnÞ : ð11Þ In our experiments we have used a2 = 1. This is not an optimal method since it requires 2 functionand 3 derivative-evaluation per cycle. Another eighth order method is the optimal scheme due to Neta and Petković [14]. It is based on Kung and Traub’s optimal fourth order method [8] and inverse interpolation. yn 1⁄4 xn un; zn 1⁄4 xn f ðynÞ f 0 n 1 1⁄21 f ðynÞ=fn 2 ; xnþ1 1⁄4 xn fn f 0 n þ cnf 2 n dnf 3 n ; ð12Þ where dn 1⁄4 1 1⁄2f ðynÞ fn 1⁄2f ðynÞ f ðznÞ yn xn f ðynÞ fn 1 f 0 n 1 1⁄2f ðynÞ f ðznÞ 1⁄2f ðznÞ fn zn xn f ðznÞ fn 1 f 0 n ; cn 1⁄4 1 f ðynÞ fn yn xn f ðynÞ fn 1 f 0 n dn1⁄2f ðynÞ fn : ð13Þ Neta’s 16th order family of optimal methods [15] is given by yn 1⁄4 xn un; zn 1⁄4 yn f ðynÞ f 0 n fn þ bf ðynÞ fn þ ðb 2Þf ðynÞ ; tn 1⁄4 xn fn f 0 n þ cnf 2 n dnf 3 n ; xnþ1 1⁄4 xn fn f 0 n þ qnf 2 n cnf 3 n þ qnf 4 n ; ð14Þ where cn and dn are given by (13) and Fig. 2. 2.682 2588 M. Scott et al. / Applied Mathematics and Computation 218 (2011) 2584–2599 qn 1⁄4 /ðtnÞ /ðznÞ FðtnÞ FðznÞ /ðynÞ /ðznÞ FðynÞ FðznÞ FðtnÞ FðynÞ ; cn 1⁄4 /ðtnÞ /ðznÞ FðtnÞ FðznÞ qnðFðtnÞ þ FðznÞÞ; qn 1⁄4 /ðtnÞ cnFðtnÞ qnF ðtnÞ ð15Þ and for dn = yn, zn, tn FðdnÞ 1⁄4 f ðdnÞ fn; /ðdnÞ 1⁄4 ðdn xnÞ FðdnÞ 1 f 4 n FðdnÞ : ð16Þ In our experiments we have used b = 2. 2. Numerical experiments We have used the above methods for 6 different polynomials. Some have real and some have complex coefficients. One example have only real roots and the rest have a combination of real and complex ones. All the roots are simple. In the first case we have taken the cubic polynomial x þ 4x 10: ð17Þ Clearly, one root is real (1.365230013) and the other two are complex conjugate. Note that the basin of attraction of each root is larger for Halley’s method than Newton’s, see Fig. 2. We have shown the results for King’s method using the parameter b 1⁄4 2. The results are not much better for several other values of b we tried. The basins of attraction for the optimal Kung–Traub method is better than any of the King’s method (notice the second quadrant of Fig. 3). Murakami’s fifth order method gives basins of attraction similar to Newton’s, see Fig. 4. On the other hand, Neta’s sixth order method which is based on King’s method and uses b = 1/2 shows some chaotic behavior in the second quadrant (Fig. 4). One sees too many points there that converge to the root in the third quadrant. This is similar to the results for King’s method. Neta–Johnson’s eighth order method has basins of attraction similar to Newton’s method, see Fig. 5. On the other hand, the optimal eighth order method (12) is much more chaotic (see Fig. 6). The 16th order optimal method (see Fig. 6) has a smaller basin of attraction for the real root, but it does not show the chaotic behavior of King’s, Neta’s sixth order and Neta–Petkovic eighth order. In general, one cannot say that increasing the order of the method will adversely affect the basins of attraction very much. In our next example, we took a quintic polynomial with real simple roots. It is clear that the best methods are Newton’s (Fig. 7), Halley’s (Fig. 7), Murakami’s (Fig. 9) and Neta–Johnson’s (Fig. 10) schemes. See Figs. 7–11. x 5x þ 4x: ð18Þ –3 –2 –1 1 2 3 –3 –2 –1 1 2 3 –3 –2 –1 1 2 3